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We use lattice Boltzmann simulations to study the ef- 
fect of shear on the phase ordering of a two-dimensional bi- 
nary fluid. The shear is imposed by generalising the lattice 
Boltzmann algorithm to include Lees-Edwards boundary con- 
ditions. We show how the interplay between the ordering 
effects of the spinodal decomposition and the disordering ten- 
dencies of the shear, which depends on the shear rate and 
the fluid viscosity, can lead to a state of dynamic equilibrium 
where domains are continually broken up and re-formed. 



PACS numbers: 83.10.Lk; 64.60.Qb; 47.n.+j 



I. INTRODUCTION 

We present numerical results for the effect of shear flow 
on the spinodal decomposition of a two-dimensional bi- 
nary fluid using lattice Boltzmann simulations. We show 
how the lattice Boltzmann algorithm can be generalised 
to allow the introduction of the Lees-Edwards bound- 
ary conditions, which are commonly used in molecular 
dynamics simulations to impose a shear flow without in- 
troducing walls. Results are presented showing how the 
competition between the ordering effects of the free en- 
ergy and the disordering effects of the shear influences the 
spinodal decomposition and piiase ordering of the fluid. 
For a recent review see Onukitl. 

When a binary fluid consisting of an equal amount of 
two components, A and B say, is rapidly cooled below the 
critical temperature it phase separates into an A-rich and 
B-rich phase. Once well-deflned domains of each phase 
are formed the typical domain size grows according to a 
power law 



R{t) - e 



(1) 



where a is the growth exponentB. a depends on the 
growth mechanism, which is dictated by the surface ten- 
sion, viscosity and diffusivity of the fluid, and the time 
elapsed after the quench. In two-dimensional systems 
diffusive Lifshitz-Slyozov growth gives a = ^ while hy- 
drodynamics can lead to faster growth with a = | . 
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The most obvious effect of shear flow on the domain 
growth is that the growing domains are elongated in the 
direction of the flow, leading to an anisotropic morphol- 
ogy. Experiments in three dimensions have shown that 
a string-like phase of thin domains oriented parallel to 
the shear can be formed in strong shearB. Such domains, 
which would normally be expected to be unstable due to 
the Rayleigh instability, appear to be stabilised by the 
shear, although very recent experiments, show that they 
can eventually break up in strong shears. 

This apparent stabilization suggests the possibility of 
a dynamic equilibrium when stretching and breaking of 
the domains as the result of the shear is balanced by 
their growth due to the thermodynamic driving force and 
to the coalescence of the domains, which can itself be 
driven by the shear. This was flrst proposed by Ohta 
and NozakiO on the basis of two-dimensional simulations 
using a cell dynamic approach. These simulations, how- 
ever, did not include hydrodynamics. 

Simulations of phase separation under shear which in- 
clude hydrodynamics are limited. Rothman performed 
early work using lattice gas cellular automata in two and 
three dimensions and was able to see the anisotropy of the 
growthB'Q. Wu et. al. undertook Langevin simulations in 
two and three dimensions and report the efzjentual forma- 
tion of a string phase in three dimensionsB. Padilla and 
Toxvaerd performed molecular dynamics simulations on 
a two-dimensional Lennard- Jones system, again jpointing 
out the anisotropic nature of the domain growthEI. In the 
simulations a peak was seen in the excess shear viscosity 
as a function of time corresponding to the increase in the 
lengths of interfaces in the system. However, there seems 
to be no evidence for a shear-induced dynamic equilib- 
rium. 

Here we simulate phase separation under shear using a 
lattice Boltzmann scheme in the same spirit as the model 
introduced by Orlandini et. al.^ which imposes phase 
separation by deflning the fljikLequilibrium as the mini- 
mum of an input free energytjllj. This method has been 
very successful in obtaping results for phase separation 
in the absence of shearll3. A particular advantage of the 
approach is that the fluid viscosity and diffusivity can be 
tuned, and this has allowed us to compare simulations for 
parameter values where diffusive or hydrodynamic phase 
separation dominates. We flnd either phases striped in 
the shear direction, or a dynamic equilibrium where the 
length scales remain approximately constant in time, de- 
pending on the relative strengths of the shear and the 
ordering. 

The lattice Boltzmann approach is described in §2. Be- 
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cause this is a lattice rather than particulate simulation 
method, it is not immediately obvious how to define Lees- 
Edwards shear boundary conditions. An approach for 
doing this is given in §3. In particular it is necessary to 
generalize the normal definition of the lattice Boltzmann 
equilibrium distribution. In §4 we define suitable mea- 
sures to characterise the anisotropic morphology of the 
spinodal decomposition patterns when shear is applied. 
The results of our simulations are contained in §5, where 
the effect of shear is compared for different fiuid viscosi- 
ties. §6 summarises the results and discusses outstanding 
questions. 



where P^js is the pressure tensor, F is a mobility parame- 
ter, ji is the chemical potential for the density difference 
and 5 is the Kronecker delta. The physical motivation 
for these constraints is twofold; firstly to ensure the cor- 
rect form of the macroscopic equations of motion and 
secondly to reproduce the correct thermodynamics of the 
binary mixture in equilibrium as discussed in more detail 
below. 

Taylor-expanding the evolution equations (||) and (^) 
to second order in the derivatives gives-Lhe macroscopic 
equations of motion for the binary fiuidt^. These are the 
continuity equation for the total density 



dtn + docTiUoc = 0, 



(9) 



II. THE LATTICE BOLTZMANN APPROACH 

The starting point for lattice Boltzmann simulationsEl 
is the evolution equation, discrete in space and time, 
for a set of distribution functions, /i, each associated 
with a velocity vector, v^. For the sake of simplicity 
we consider a pkigle relaxation time, the so-called BGK 
approximationllj. The evolution equation for the {fi} is 



/i(x + ViAt,t + At) 



/,(x,t) = — (/f-/,), (2) 
n 

where x is a lattice point. At is the time step, and ViAt 
is normally constrained to be a lattice vector. The relax- 
ation time is ri and is the equilibrium distribution. 
For a two-component system a second, equivalent equa- 
tion is also needed 

^,(x + v,At, t + At) - ^,(x, t) = —{g'l - Qi). (3) 

Physical quantities are defined as moments of the dis- 
tribution functions. To model the isothermal fiow of a 
binary mixture of components A and B, we choose 



nu, 



^9i = ^^ 



(4) 



where n is the total density field, u is the velocity field 
and (f is the field corresponding to the difference in the 
density of components A and B. 

We require mass conservation for both components and 
momentum conservation for the bulk. This is equivalent 
to constraining the equilibrium distributions to obey 



nu. 



(5) 



We also need to define higher-order moments of the 
equilibrium densities. The choice for these moments is 
withi-Qpthe free energy lattice Boltzmann scheme used 
herdiaM 
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(6) 

(7) 
(8) 



a convection-diffusion equation governing the evolution 
of the density difference 

dt^ + dai^U^) = UJ2 (fVV - S/3 [^do^Pa^y) , (10) 

and, in the incompressible limit, the incompressible 
Navier Stokes equations for a non-ideal system 

ndtUa + nUf3df3Ua = -df3Paf3 + V^ii^ + 0{d^) (11) 

where cji,2 = ti,2 — At/ 2 and the viscosity is given by 
u = ncoi/S. 

The thermodynamic fields entering the simulation are 
the pressure tensor and the chemical potential which fol- 
low from the free energy of the system. We consider the 
free energy of a simple binary fiuid. A-A and B-B inter- 
actions are zero, but there is an A-B repulsion Xtiatib 
where ua and ub are the number densities of A- and 
B-particles, respectively, and A is a parameter describing 
the interaction strength. This system can be described 
by the Landau free energy functional 

* = j dv[i,{^,n,T) + ^ (y^f } (12) 

where T is the temperature and z^: is a measure of the 
excess interface free energy (surface tensioni— The free 
energy density of the homogeneous system iaij 



(13) 



+ -{n + ^) In l—^ 

fn-ip 



For temperatures greater than a critical temperature 
Tc = A/ 2 the system remains in a single phase. For 
T < Tc there is phase separation into two states with 

From the free energy (|l^) we derive the local chemi- 
cal potential /i as the functional derivative of the total 
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free energy ^ with respect to the concentration difference 
field (^(x) 

= ^—r^ = + - In - nV^LD. (14) 

Equilibrium corresponds to n, T) = 0. 

The derivation of the pressure tensor |k slightly more 
involved and is discussed in Appendix At3. We obtain 

= {nT ^ (p/J.^{x))Sal3 

-\-K,{da^di3(p — -d^(pd^(pSa/3 — (fd^d^if) 
= {nT ^ (pfl{x))Sal3 

^K{da(pdf3(p - ^d^cpd^ipSa/s) (15) 

where the first term is the ideal gas pressure, the sec- 
ond term is the osmotic pressure with ja^ = dcjyt/j and 
the third term is related to the surface tension. The os- 
motic prepjpe was omitted in the original definition of 
the modeOO. The chemical potential and pressure ten- 
sor are input to the lattice Boltzmann scheme through 
equations (||) and (^). In equilibrium the simulated fluid 
minimises the free energy (|l2|). 

It remains only to define the equilibrium distributions 
and introduced in the evolution equations (||) and 
(H) . Normally an expansion to second order in the veloc.- 
ities is sufficient to reproduce the constraints (||) - (^)li2l. 
However, this ceases to be the case when Lees-Edwards 
shear boundary conditions are introduced. In the next 
section we discuss how the equilibrium distribution can 
be defined to allow the use of Lees-Edwards boundary 
conditions. 



III. SHEAR BOUNDARY CONDITIONS 

Possibly the easiest way to introduce shear flow in a 
lattice Boltzmann simulation is to include walls moving 
in a lattice direction. Even for a wall with neutral wet- 
ting, however, phase separation is strongly enhanced at 
the walls and the wall effects easily dominate the phase 
separation process for all but the largest systems. The 
effect of walls on phase separation is an interesting phe- 
nomenon in its own right, but it is not the process we are 
interested in studying here. 

The problem caused by explicit walls can be overcome 
in a relatively simple and efficient manner by introducing 
a Klein-bottle symmetry to the lattice. This is done by 
forcing the fluid to have a given velocity along one line 
in the direction of the shear flow. In a one-component 
mixture this induces a linear velocity profile. For a two- 
component mixture, however, the dynamics are influ- 
enced by the V-shaped velocity profile at the forcing line 



because of the non-local interactions. We used this al- 
gorithm to produce preliminary results but it has no ad- 
vantages over the method derived below. 

A more regular shear flow can be produced by ex- 
tending the idea of Lees-Edwards bpojudary conditions, 
widely used in Molecular Dynamicalj, to lattice Boltz- 
mann simulations. Briefly, Lees and Edwards simulated 
shear boundary conditions for a shear in the x-direction 
in a simulation box of dimensions (X, Y) by introducing 
periodic boundary conditions in the ^/-direction. Parti- 
cles that left the box at the lower boundary for position 
{x^y = 0) reappeared at the upper boundary at position 
{x-\-ut{ mod X), y = Y) with a velocity that was changed 
hy V ^ V -\- u. 

To implement this idea for lattice Boltzmann simula- 
tions we are faced with two difficulties. Firstly the densi- 
ties are defined on a lattice and the Lees-Edwards bound- 
ary conditions lead to densities defined between the lat- 
tice points. Secondly we need to impose a Galilean trans- 
formation for the densities which are streamed across the 
lattice. 

The non-fitting of the lattice is relevant for both the 
streaming and for the calculation of derivatives at y = 1 
and y = Y. We solve this problem by a linear interpola- 
tion scheme. For any density we define 

f[x, y = 0] = {l- R{ut))f[x + I{ut), y = Y] 

^R{ut)f[x + I{ut) + 1, ^ = r] (16) 

where I{z) is the largest integer with I{z) < z and 
R{z) = z — I{z). If we pass the break in the lattice 
from the other side we define similarly 

f[x, y = y + 1] = (1 - R{ut))f[x - I{ut), y = l] 

+R{ut)f[x-I{ut)-l,y=l]. (17) 

These formulae are used both for the streaming of the 
Galilean-transformed Boltzmann densities, /i, and for 
the calculation of density gradients. 




FIG. 1. Our numbering of the velocity vectors in a 
nine- velocity model. 

It is rather more difficult to see how the Galilean trans- 
formation should be defined. Let us consider the special 
case of a two-dimensional, nine-velocity model where the 
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velocities are numbered as indicated in Figure |l|. We need 
to perform a Galilean transformation on the {5,2,6} and 
the {7,4,8} velocities as these will carry mass and mo- 
mentum across the boundaries. To define the transfor- 
mation we demand mass and ^-momentum conservation 



Hp = /5 + /2 + /6 = /s + /2 + /e 

an appropriate change in the x-momentum 

ih - k) - (/5 - /e) = 
and conservation of the local pressure 

Pxx = -^(^^^x - npUa:f 
Tip 



(18) 



(19) 



-^fei-rip - UpUxf) 
1 

Up 



(20) 



where the prime denotes the transformed quantities. 

This system of equations can be solved to give a unique 
solution for the Galilean-transformed densities f[ 



/2 = /2 + 2(/5 - fQ)u - UpU^, 

3 11 

/s = /5 + (- - + i^h)"^ - 

/6 = /6 + (-U + U + ?/6)^- 



(21) 
(22) 

(23) 



This definition can be extended to a Galilean transfor- 
mation for all densities and, equivalently, to a transfor- 
mation in different lattice directions. 

In order for this transformation to make sense we need 
to make sure that equation (|l^) is consistent with the def- 
inition of the equilibrium distribution, in Eqn. (H) i.e. 
that an equilibrium distribution for a velocity u Galilean 
transformed by a velocity Aix is equal to the equilibrium 
distribution for velocity u -|- Atx. It is conventional to 
define the equilibrium distribution as a polynomial in 
second order in u. A generic expansion is 

= AcrTl + BcrTlUaVic, + CaUU^ 

-^D^nUaUpVic^Vip + Gaa(3nVic,Vif3 (24) 

where A^r, D(j^ GaafB are constants that have the 

absolute value of the corresponding velocity vector a = 
\vi\ as an index. However, substituting (24) into (|l9| ) 
shows that this equation is not satisfied in equilibrium. 
In practice this leads to a step in the Ux profile at the 
boundary. 

There is, however, no a priori reason to use a second- 
order expansion in the velocity for the equilibrium dis- 
tribution. All that is needed for a valid equilibrium dis- 
tribution is that (||)~(B hold and that the distribution 
obeys the conditions (|2l])-(|2^). 



Let Tap = Pafd/n. Then, if we require. 



Ji Is 

f? + fS + fl 



(25) 



/5 + /e - if 5 + /e + f2)iT.. + u,u,)=0, (26) 

/5° + /S - (/S + /S + f?)iTyy + %%) = 0, (27) 

fs + fr - ifs + fi + f!){T.. + u,u,) = 0, (28) 
/e + fr - (/e + /r + fsKTyy + UyUy) = 0. (29) 

29|) , together with (|5|)-(||) are a completely deter- 
mined set of equations with the solution 



f! 

fi 
fi 
fi 
fi 



^yy 



-n(TxX + + U^){1 - Tyy - Uy), 

n(Tyy^Uy^ul)(l-TxX-ul), 

n(Txx -Ux^ ul){l - Tyy - ul), 



2 

1 

2 

\<Tyy -Uy^ U^) (l " T^^ " ) 

~^^(Txy -r j-xx-^yy ^ -^yy 
,2 



" ^xx^yy Pyyi'^x ^x) 



^TxxiUy + Uy) + UxUyil ^Ux^Uy^ UxUy)) , 
'^( '^Xy H~ TxxTyy -\- Tyyiy Ux H~ x^ 



^ ■ - \ -^y ' - 0.0.- yy I -^yy 

-^Txx{Uy + Uy) - UxUy{l -Ux^Uy- UxUy)) 

-^^(Txy -r j-xx-^yy ~r -^yy 
,2 



' TxxTyy + Tyy( Ux Ux) 



U. 



Uy + Uy) + UxUy{l 

~^^{ ^xy ^xx^yy ^yyi^x H~ U^) 



y 



- u. 



Uy))^ 



^Txx{-Uy + U^) - UxUy{l -^Ux-Uy - UxUy)). 

For this equilibrium distribution 

f _ f rp 

J 5 J6 _ ^, J -'-xy 



fi + fl + fl 



(30) 



"-yy 



which is consistent with the Galilean transformation (|T^). 

For a two-component system we similarly define the 
using 



9'2+9'o+9i 



and imposing 



(31) 



(32) 



where ^ is a free parameter that can be used to improve 
stability (we choose £ = 1). Solving equations (|3l] ) and 
and dq)-(§) gives 
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5? = l/2{(^- 

<?5° = l/4{(2- 
<?6° = l/4{(2- 
l/4{(2- 
l/4{(2- 



5? 
9% 



Ux)Vii + (1 + - u\)i^Ux}, 

- Uy)T^ - {1 - Uy - ul)ipUy}, 

- Ux + Uy)r/1 + (1 + iia, + Uy)(pUxUy}, 
Ux + Uy)Tll - (l-Ux^ Uy)(pUxUy}, 
Ux - Uy)Tll {l-Ux- Uy)(fUxUy}, 

- Ux - Uy)T/l - {l^Ux - Uy)(pUxUy}. 



The macroscopic flow equations are unaffected by the 
choice of the further constraints (p5|)-(29) and (^-(|3^) 
or by the detailed structure of the equilibrium distri- 
butions. Therefore, these alterations in the model can 
change the numerical stability and the behaviour of quan- 
tities like the spurious velocities, but they leave the evo- 
lution of the macroscopic quantities unaffected, at least 
to second order in the derivatives. 



IV. MEASURES FOR NON-ISOTROPIC 
PATTERNS 

To characterise the features of phase separation under 
shear it is necessary to construct measures for the length 
scales of the sheared systems which will in general be 
anisotropic. Measures that are based on Fourier trans- 
forms cannot be easily used for sheared systems because 
the system is no longer periodic. 

Length scales derived from derivatives do not require 
periodicity. Derivatives need to be evaluated for the al- 
gorithm and are readily available. We define a tensor 



daf3 



(33) 



where is the symmetric discrete derivative in direction 
a. Because the tensor is symmetric it can be diagonalised 
to give two eigenvalues Ai, A2 and an angle 9^ 



Ai = 



■V dxx H~ dyy 

M = z 



r = tan"^ 



{dx 



- d 12 

4 ^xy> 



(dxx dyy) 



■ (P 



dxy — A2 



(34) 

(35) 
(36) 



The two eigenvalues give two orthogonal length scales 
1 1 



\2{t)Lw ' 



(37) 



where Ly^ is the interface width. It appejars because d^^ 
scales inversely with the interface widthllj. L^, used as 
a constant here, could in principle be anisotropic. That 



this anisotropy is not a strong effect can be seen by com- 
paring these length scales with scales that are explicitly 
independent of the interface width. 

One such measure is related to the lengths of the in- 
terfaces in the system. The interface can be represented 
by a set of contours. These contours consist of small line 
segments U and the length of the interface can be written 



(38) 



In order to extract the preferred direction of the interface 
we define the vector 



D = R- 



The operator R is defined by 



R{x) = \x\ 



cos(2i9) \ 
sin(2i9) J 



(39) 



(40) 



where is the angle between the argument of R and the 
X-axis. 

5 is a vector that is zero for isotropic closed contours 
and which points in the average direction of the inter- 
face for non-isotropic closed contours. Two length scales 
and an angle that correspond to the intuitive result for 
oriented rectangular objects can be defined from these 
measures 



XY 



Li^\D\ 



= cos 



Rn 



x.D 



XY 



L, 



\D\ 



(41) 



(42) 



Thus we have defined two independent sets of measures 
for the structure of non-isotropic patterns that will now 
be used to examine spinodal decomposition under shear. 



V. SIMULATION RESULTS 

For all the simulations we used a total density n = 2, 
an interaction parameter A = 1.1, which corresponds to 
a critical temperature Tc = 0.55, and a temperature T = 
0.5. The equilibrium values of the order parameter were 
then (po = ±1. The mobility was F = 2, the relaxation 
time for the order parameter in Eqn. (||) was r2 = 1 and 
the interface free energy parameter was k = 0.002, which 
corresponds to an interface width of approximately three 
lattice spacings. The relaxation parameter for the total 
density Eqn. (|^), ri, was varied: ri = 100 gave a high 
viscosity and ri = 1 an intermediate viscosity. 

The shear transformation, S, is defined as 



5 



y 



(43) 



Shear flow applied to a system undergoing spinodal de- 
composition stretches the original pattern. This effect is 
only relevant once the deformation caused by the flow is 
of the same order or larger than the deformation caused 
by the coarsening process. This requires 



jt > 1. 



(44) 



We therefore expect to observe the effect of the shear flow 
for t > 1/7. 




time=0 



time=l 




time=2 



time=3 




time=4 



time=5 




FIG. 2. Applying shear (with shear rate 7 = 1) to ^ 
without internal dynamics leads to homogenization. 



, system 



To help understand the effect of shear-flow on a phase- 
separating system let us first consider a pattern without 
any internal dynamics that undergoes a shear transfor- 
mation. This transformation is illustrated in Figure 0, 



where we start from a frozen spinodal decomposition pat- 
tern and show successive iterations of a shear transfor- 
mation with 7 = 1. 



nme=Z3U 




10000 20000 



time time 
(b) (c) 
FIG. 3. (a) Spinodal decomposition under shear for a high 
viscosity binary fluid (n = 100, L^: = 256, = 128). The 
high viscosity suppresses internal hydrodynamic degrees of 
freedom. The shear rate is 7 = 0.004 which corresponds to 
a shear time ts = 250. (b) Variation of the orientation (in 
degrees) of the pattern with time, (c) Variation of the length 
scales (in arbirtary units) with time. 

The structure develops an orientation that slowly 
aligns with the shear direction while the stretching in- 
creases the length of the domains along the shear. Once 



6 



the width of the domains is smaher than the original 
width of the interface the system is effectively a homoge- 
neous mixture. 

This effect is known as shear-induced mixing. It can be 
observed in the lattice Boltzmann fluids if the stretching 
effect of the shear flow is much faster than the growth 
of the domains via diffusion or flow. Numerically this 
can be achieved by choosing a very low mobility and a 
high viscosity. Phase separation is suppressed because of 
the mixing properties of the shear flow unless the phase- 
separating structure is aligned with the shear direction. 
For finite lattices we sometimes observe at much later 
times a nucleation of complete stripes that span the sys- 
tem and are periodic in the shear direction. The time 
required to form these stripes depends on the system size 
and it seems reasonable to assume that this phenomenon 
does not occur in infinite systems. 

We now consider a high viscosity fluid (n = 100) in 
which diffusive but not hydrodynamic modes are impor- 
tant. The internal dynamics leads to domain coarsening 
and can also prevent a complete mixing of the system. 
Figure ^ shows the spinodal decomposition pattern of 
the high- viscosity binary mixture. For very short times 
{t < 300 ~ 7~^) we observe the familiar spinodal decom- 
position pattern. It is, however, coarsening in a new way 
via shear flow-induced collisions of the domains. This 
process enhances domains oriented in the collision direc- 
tion. Then for 300 < t < 1000 the flow slowly turns 
the striped pattern and stretches it. At t ^ 1000 the 
rupturing of domains starts to be important and for 
1000 < t < 15000 there is a continuous stretching and 
rupturing that effectively stops the phase ordering pro- 
cess. At t ^ 15000 the system developes stripes that 
span the system. Because periodic stripes are unaffected 
by the shear flow if they are completely aligned with it 
the system now grows via the diffusion mechanism. 

This evolution can be followed more quantitatively by 
measuring the orientation angle and the length scales de- 
fined in Section IV. Figure ^3 shows the angle of orienta- 
tion to the X-axis measured by 0^ (Eqn. ^6|) and ^° (Eqn. 
^). The two different measures for the angle agree very 
well. The pattern tilts at very early times {t < 2000) and 
then slowly aligns with the direction of the shear flow as 
periodic stripes are created. 

The graph in Figure ||c shows the length-scales Ri 2 
defined in Eqns. ( |37| ) and the length scales 2 defined 
m Eqns. (|ll). We very clearly see a separation of length 
scales and a good agreement of the two different mea- 
sures. A minimum of the larger length scale at t ~ 17000 
indicates the creation of periodic stripes spanning the 
system. After this time the growth of domains is no 
longer hindered by the continual breaking of stretched 
domains. 

We now turn to consider a system with a lower viscos- 
ity that allows for a hydrodynamic response of the do- 
mains to the shear flow. Results are presented in Figure 
^. It is immediately obvious that the pattern differs 




time 



(b) (c) 
FIG. 4. (a) Spinodal decomposition under shear for a 
medium viscosity binary fluid (ti = l,Lx = 256, = 128). 
The effect of internal flow causes the domains to remain at 
an angle to the shear direction. The shear rate is 7 — 0.004 
which corresponds to a shear time ts = 250. (b) Variation 
of the orientation (in degrees) of the pattern with time, (c) 
Variation of the length scales (in arbitrary units) with time. 

from that in Figure |[ The final state does not sim- 
ply consist of periodic stripes, but of dynamic structures 
that are constantly stretched, broken and deformed by 
the flow. At least on this time scale a state of dynamic 
equilibrium is reached where the ordering effects of the 
spinodal decomposition balance the disordering effects of 
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the shear. 

The quantitative measures in Figures ^ and ||c show 
that after initial fluctuations the orientation of the pat- 
tern converges to a value that fluctuates about a finite 
angle to the shear direction. This phenomenon is similar 
to the behaviour of a singk. sheared drop that lies at a 
finite angle to a shear flowll3. The graph of length scales 
again shows a very clear distinction between the large and 
small length scales. Strong oscillations are seen. These 
may be finite size effects because the system is so small 
and contains only a few domains. However such oscil- 
lations-bave been seen in experimentsu and ina model 
systeniE2l. 




time 



(b) (c) 



FIG. 5. (a) Spinodal decomposition under shear for a high 
viscosity binary fluid (n = 100, — 512, Ly = 512). The 
high viscosity suppresses internal hydrodynamic degrees of 
freedom. The shear rate is 7 = 0.0001 which corresponds to 
a shear time ts = 1000. (b) Variation of the orientation (in 
degrees) of the pattern with time, (c) Variation of the length 
scales (in arbitrary units) with time. 

We have, so far, considered strong shear flow. Let us 
now consider the same viscosity, where both diffusive and 
hydrodynamic flow is possible, but lower the shear rate 
so that the early time spinodal decomposition is unaf- 
fected by the flow. In Figure || the spinodal decomposi- 
tion for a shear rate 7 = 0.0001 is shown for a system 
with Ti = 1. For times t < I/7 = 10000 we see the typ- 
ical spinodal decomposition pattern for these viscosities. 
Hydrodynamic flow leads to circular domains which then 
grow through the slower diffusive mechanism. After this 
time, the stretching of the domains dominates over the 
domain growth and the pattern becomes non-isotropic. 
By t ^ 10000 the pattern comprises large-stripe like do- 
mains together with the nested pattern of drops within 
drops in the large domains. As the large domains are 
stretched, the drops inside them coalesce with the walls 
and slowly the stripes are cleaned of the small included 
drops. 

These results also clearly show up in the measurements 
given in Figure ||. After t > 10000 the orientation slowly 
converges towards a tilting angle 6> ~ 7°, the long and 
short length scales split and the ^ R° ^ ^ 
growth law breaks down. In the R^ measure derived 
from the number of domains we see a slight increase from 
the normal growth law corresponding to the process of 
shear cleaning the stripes from drops. 



VI. CONCLUSIONS 

In this paper we have investigated the effects of shear 
flow on systems undergoing spinodal decomposition. In 
order to study these systems we introduced an extension 
to the lattice Boltzmann algorithm that allows simula- 
tion of shear flow problems with Lees-Edwards boundary 
conditions. We find that the effect of shear flow on spin- 
odal decomposition depends strongly on the viscosity of 
the fluid. Systems with a very high viscosity tend to or- 
der in the shear direction, whereas systems with a lower 
viscosity arrive at a dynamic stationary state where the 
domains lie at a finite angle to the shear direction. 

One of the problems in simulating spinodal decompo- 
sition under shear is that the shear flow induces long- 
range correlations much faster than for un-sheared sys- 
tems so that larger lattice sizes are required to exam- 
ine long-time behaviour. Therefore there remain many 
unexplored problems concerning the structure of spin- 
odal decomposition under shear. For example, it would 
be interesting to investigate the transition between the 
sheared and non-sheared patterns for different viscosities 
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and to ask whether the late-time decomposition patterns 
are statisticahy independent of an initial shear. 



APPENDIX A 

We show how the full pressure tensor (|l5|) is derived. 
The pressure of a homogeneous system is defined as the 
volume derivative of the free energy. Writing the full 
volume dependence of the densities n = N/V and (p = 
{Na — Nb)/V explicitly we see that: 



N Na-Nb \ 
V J 
« A/ / Na-Nb 



(45) 



For a non-homogeneous system the pressure is no 
longer a scalar but a tensor. The correct form of the pres- 
sure tensor can be derived from a Lagrangian expression 
for the free energy which is minimized in equilibrium 



^ - {Na - Nb)) ^ [ n-N). (46) 
Jv Jv 

To obtain differential equations for the equilibrium we 
evaluate the Euler-Lagrange equations and get 



(47) 
(48) 



We multiply these equations with 9/3 and 9/3 n, respec- 
tively and sum the equations. Remembering that ji^ and 
fin are constants, this yields 



df3{^l^^ + n/in) = -9c(#a/3 



-i^{da(pdf3(p - -d^cpd^cpSafd)}- (49) 



We then substitute the expressions for the chemical po- 
tentials (47) and (48) into ( p9| ) and subtract the right- 
hand side from the left-hand side to derive a tensor a 
that has a zero divergence 

^K{da(pd(3(p - ^d^(pd^(pSai3 - (fd^d^ipda/s)) - (50) 

For a uniform system = P5a(3 reduces to the homo- 
geneous pressure. The divergence of the pressure tensor 
must vanish in equilibrium. We therefore identify Gcx(3 
with the pressure tensor Pq;/3. 
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